mini_starmap

# to automatically save figures in save_dir set save_plot to TRUE
temp_dir = '../results/' #'~/Temp/'
myinstructions = createGiottoInstructions(save_dir = temp_dir,
                                          save_plot = FALSE, 
                                          show_plot = F)
#> 
#>  no external python path or giotto environment was specified, will check if a default python path is available 
#> 
#>  A default python path was found:  /opt/conda/bin/python3  and will be used
#> 
#>  If this is not the correct python path, either
#> 
#>  1. use installGiottoEnvironment() to install a local miniconda python environment along with required modules 
#> 
#>  2. provide an existing python path to python_path to use your own python path which has all modules installed

# giotto object 
expr_path = system.file("extdata", "starmap_expr.txt.gz", package = 'Giotto')
loc_path = system.file("extdata", "starmap_cell_loc.txt", package = 'Giotto')
starmap_mini <- createGiottoObject(raw_exprs = expr_path,
                                   spatial_locs = loc_path,
                                   instructions = myinstructions)
#> Consider to install these (optional) packages to run all possible Giotto commands for spatial analyses:  SPARK multinet RTriangle FactoMiner
#>  Giotto does not automatically install all these packages as they are not absolutely required and this reduces the number of dependencies

# show instructions associated with giotto object (starmap_mini)
showGiottoInstructions(starmap_mini)
#> $python_path
#> [1] "/opt/conda/bin/python3"
#> 
#> $show_plot
#> [1] FALSE
#> 
#> $return_plot
#> [1] TRUE
#> 
#> $save_plot
#> [1] FALSE
#> 
#> $save_dir
#> [1] "../results/"
#> 
#> $plot_format
#> [1] "png"
#> 
#> $dpi
#> [1] 300
#> 
#> $units
#> [1] "in"
#> 
#> $height
#> [1] 9
#> 
#> $width
#> [1] 9
#> 
#> $is_docker
#> [1] FALSE

filterDistributions(starmap_mini, detection = 'genes')

filterDistributions(starmap_mini, detection = 'cells')

filterCombinations(starmap_mini,
                   expression_thresholds = c(1),
                   gene_det_in_min_cells = c(50, 100, 200),
                   min_det_genes_per_cell = c(20, 28, 28))

starmap_mini <- filterGiotto(gobject = starmap_mini,
                           expression_threshold = 1,
                           gene_det_in_min_cells = 50,
                           min_det_genes_per_cell = 20,
                           expression_values = c('raw'),
                           verbose = T)
#> Number of cells removed:  95  out of  4000 
#> Number of genes removed:  0  out of  28
starmap_mini <- normalizeGiotto(gobject = starmap_mini, scalefactor = 6000, verbose = T)
#> 
#>  first scale genes and then cells
starmap_mini <- addStatistics(gobject = starmap_mini)

starmap_mini <- runPCA(gobject = starmap_mini, method = 'factominer')
#> hvg  was not found in the gene metadata information, all genes will be used
#> Warning in runPCA_factominer(x = t_giotto(expr_values), scale = scale_unit, :
#> ncp > ncol(x), will be set to ncol(x)
screePlot(starmap_mini, ncp = 30)
#> PCA with name:  pca  already exists and will be used for the screeplot

plotPCA(gobject = starmap_mini)

# 2D umap
starmap_mini <- runUMAP(starmap_mini, dimensions_to_use = 1:8)
plotUMAP(gobject = starmap_mini)

# 3D umap
starmap_mini <- runUMAP(starmap_mini, dimensions_to_use = 1:8, name = '3D_umap', n_components = 3)
plotUMAP_3D(gobject = starmap_mini, dim_reduction_name = '3D_umap')
#> create 3D plot

# 2D tsne
starmap_mini <- runtSNE(starmap_mini, dimensions_to_use = 1:8)
plotTSNE(gobject = starmap_mini)

starmap_mini <- createNearestNetwork(gobject = starmap_mini, dimensions_to_use = 1:8, k = 25)
starmap_mini <- doLeidenCluster(gobject = starmap_mini, resolution = 0.5, n_iterations = 1000)

# 2D umap
plotUMAP(gobject = starmap_mini, cell_color = 'leiden_clus', show_NN_network = T, point_size = 2.5)

# 3D umap
plotUMAP_3D(gobject = starmap_mini, dim_reduction_name = '3D_umap', cell_color = 'leiden_clus')
#> create 3D plot
#>  center label is not clear to see in 3D plot
#>  You can shut it down with show_center_label = F

# 2D umap + coordinates
spatDimPlot(gobject = starmap_mini, cell_color = 'leiden_clus',
            dim_point_size = 2, spat_point_size = 2.5)

# 3D umap + coordinates
spatDimPlot3D(gobject = starmap_mini, cell_color = 'leiden_clus', dim_reduction_name = '3D_umap')
#>  center label is not clear to see in 3D plot
#>  You can shut it down with show_center_label = F
#> Warning: 'layout' objects don't have these attributes: 'NA'
#> Valid attributes include:
#> 'font', 'title', 'uniformtext', 'autosize', 'width', 'height', 'margin', 'paper_bgcolor', 'plot_bgcolor', 'separators', 'hidesources', 'showlegend', 'colorway', 'datarevision', 'uirevision', 'editrevision', 'selectionrevision', 'template', 'modebar', 'meta', 'transition', '_deprecated', 'clickmode', 'dragmode', 'hovermode', 'hoverdistance', 'spikedistance', 'hoverlabel', 'selectdirection', 'grid', 'calendar', 'xaxis', 'yaxis', 'ternary', 'scene', 'geo', 'mapbox', 'polar', 'radialaxis', 'angularaxis', 'direction', 'orientation', 'editType', 'legend', 'annotations', 'shapes', 'images', 'updatemenus', 'sliders', 'colorscale', 'coloraxis', 'metasrc', 'barmode', 'bargap', 'mapType'

# heatmap and dendrogram
showClusterHeatmap(gobject = starmap_mini, cluster_column = 'leiden_clus')

showClusterDendrogram(starmap_mini, h = 0.5, rotate = T, cluster_column = 'leiden_clus')

gini_markers = findMarkers_one_vs_all(gobject = starmap_mini,
                                                  method = 'gini',
                                                  expression_values = 'normalized',
                                                  cluster_column = 'leiden_clus',
                                                  min_genes = 20,
                                                  min_expr_gini_score = 0.5,
                                                  min_det_gini_score = 0.5)
#> 
#>  start with cluster  1 
#> 
#>  start with cluster  2 
#> 
#>  start with cluster  3 
#> 
#>  start with cluster  4 
#> 
#>  start with cluster  5 
#> 
#>  start with cluster  6 
#> 
#>  start with cluster  7 
#> 
#>  start with cluster  8

# get top 2 genes per cluster and visualize with violinplot
topgenes_gini = gini_markers[, head(.SD, 2), by = 'cluster']
violinPlot(starmap_mini, genes = topgenes_gini$genes, cluster_column = 'leiden_clus')
#> These genes have duplicates: 
#>  Egr1 Vip Vip Vip Pvalb Vip Sst Sst Egr2

# get top 6 genes per cluster and visualize with heatmap
topgenes_gini2 = gini_markers[, head(.SD, 6), by = 'cluster']
plotMetaDataHeatmap(starmap_mini, selected_genes = topgenes_gini2$genes,
                    metadata_cols = c('leiden_clus'))

clusters_cell_types = c('cell A', 'cell B', 'cell C', 'cell D',
                        'cell E', 'cell F', 'cell G', 'cell H')
names(clusters_cell_types) = 1:8
starmap_mini = annotateGiotto(gobject = starmap_mini, 
                              annotation_vector = clusters_cell_types, 
                              cluster_column = 'leiden_clus', 
                              name = 'cell_types')
# check new cell metadata
pDataDT(starmap_mini)

# visualize annotations
spatDimPlot(gobject = starmap_mini, cell_color = 'cell_types', 
            spat_point_size = 2, dim_point_size = 2)

dimGenePlot3D(starmap_mini,
              dim_reduction_name = '3D_umap',
              expression_values = 'scaled',
              genes = "Pcp4",
              genes_high_color = 'red', genes_mid_color = 'white', genes_low_color = 'darkblue')

spatGenePlot3D(starmap_mini,
               expression_values = 'scaled',
               genes = "Pcp4",
               show_other_cells = F,
               genes_high_color = 'red', genes_mid_color = 'white', genes_low_color = 'darkblue')

starmap_mini <- createSpatialGrid(gobject = starmap_mini,
                              sdimx_stepsize = 200,
                              sdimy_stepsize = 200,
                              sdimz_stepsize = 20,
                              minimum_padding = 10)
showGrids(starmap_mini)
#> The following grids are available:  spatial_grid
#> [1] "spatial_grid"

# visualize grid
spatPlot2D(gobject = starmap_mini, show_grid = T, point_size = 1.5)

plotStatDelaunayNetwork(gobject = starmap_mini, maximum_distance = 200, 
                        method = 'delaunayn_geometry')
#> Warning: Ignoring unknown parameters: binwidth, bins, pad

starmap_mini = createSpatialNetwork(gobject = starmap_mini, minimum_k = 2, 
                                    maximum_distance_delaunay = 200, 
                                    method = 'Delaunay', 
                                    delaunay_method = 'delaunayn_geometry')
starmap_mini = createSpatialNetwork(gobject = starmap_mini, minimum_k = 2, 
                                    method = 'kNN', k = 10)
showNetworks(starmap_mini)
#> The following images are available:  Delaunay_network kNN_network
#> [1] "Delaunay_network" "kNN_network"

# visualize the two different spatial networks  
spatPlot(gobject = starmap_mini, show_network = T,
         network_color = 'blue', spatial_network_name = 'Delaunay_network',
         point_size = 2.5, cell_color = 'leiden_clus')

spatPlot(gobject = starmap_mini, show_network = T,
         network_color = 'blue', spatial_network_name = 'kNN_network',
         point_size = 2.5, cell_color = 'leiden_clus')

km_spatialgenes = binSpect(starmap_mini)
#> 
#>  This is the single parameter version of binSpect
#>  1. matrix binarization complete 
#> 
#>  2. spatial enrichment test completed 
#> 
#>  3. (optional) average expression of high expressing cells calculated 
#> 
#>  4. (optional) number of high expressing cells calculated
spatGenePlot(starmap_mini, expression_values = 'scaled', 
             genes = km_spatialgenes[1:4]$genes,
             point_shape = 'border', point_border_stroke = 0.1,
             show_network = F, network_color = 'lightgrey', point_size = 2.5,
             cow_n_col = 2)

rank_spatialgenes = binSpect(starmap_mini, bin_method = 'rank')
#> 
#>  This is the single parameter version of binSpect
#>  1. matrix binarization complete 
#> 
#>  2. spatial enrichment test completed 
#> 
#>  3. (optional) average expression of high expressing cells calculated 
#> 
#>  4. (optional) number of high expressing cells calculated
spatGenePlot(starmap_mini, expression_values = 'scaled', 
             genes = rank_spatialgenes[1:4]$genes,
             point_shape = 'border', point_border_stroke = 0.1,
             show_network = F, network_color = 'lightgrey', point_size = 2.5,
             cow_n_col = 2)

silh_spatialgenes = silhouetteRank(gobject = starmap_mini) # TODO: suppress print output
spatGenePlot(starmap_mini, expression_values = 'scaled', 
             genes = silh_spatialgenes[1:4]$genes,
             point_shape = 'border', point_border_stroke = 0.1,
             show_network = F, network_color = 'lightgrey', point_size = 2.5,
             cow_n_col = 2)

# 1. calculate spatial correlation scores 
ext_spatial_genes = km_spatialgenes[1:20]$genes
spat_cor_netw_DT = detectSpatialCorGenes(starmap_mini,
                                         method = 'network', 
                                         spatial_network_name = 'Delaunay_network',
                                         subset_genes = ext_spatial_genes)

# 2. cluster correlation scores
spat_cor_netw_DT = clusterSpatialCorGenes(spat_cor_netw_DT, 
                                          name = 'spat_netw_clus', k = 6)
heatmSpatialCorGenes(starmap_mini, spatCorObject = spat_cor_netw_DT, 
                     use_clus_name = 'spat_netw_clus')

netw_ranks = rankSpatialCorGroups(starmap_mini, 
                                  spatCorObject = spat_cor_netw_DT, 
                                  use_clus_name = 'spat_netw_clus')
top_netw_spat_cluster = showSpatialCorGenes(spat_cor_netw_DT, 
                                            use_clus_name = 'spat_netw_clus',
                                            selected_clusters = 6, 
                                            show_top_genes = 1)

cluster_genes_DT = showSpatialCorGenes(spat_cor_netw_DT, 
                                       use_clus_name = 'spat_netw_clus',
                                       show_top_genes = 1)
cluster_genes = cluster_genes_DT$clus; names(cluster_genes) = cluster_genes_DT$gene_ID


starmap_mini = createMetagenes(starmap_mini, gene_clusters = cluster_genes, name = 'cluster_metagene')
spatCellPlot(starmap_mini,
             spat_enr_names = 'cluster_metagene',
             cell_annotation_values = netw_ranks$clusters,
             point_size = 1.5, cow_n_col = 3)

hmrf_folder = paste0(temp_dir,'/','11_HMRF/')
if(!file.exists(hmrf_folder)) dir.create(hmrf_folder, recursive = T)

# perform hmrf
my_spatial_genes = km_spatialgenes[1:20]$genes
HMRF_spatial_genes = doHMRF(gobject = starmap_mini,
                            expression_values = 'scaled',
                            spatial_genes = my_spatial_genes,
                            spatial_network_name = 'Delaunay_network',
                            k = 6,
                            betas = c(10,2,2),
                            output_folder = paste0(hmrf_folder, '/', 'Spatial_genes/SG_top20_k6_scaled'))
#> 
#>  expression_matrix.txt already exists at this location, will be overwritten 
#> 
#>  spatial_genes.txt already exists at this location, will be overwritten 
#> 
#>  spatial_network.txt already exists at this location, will be overwritten 
#> 
#>  spatial_cell_locations.txt already exists at this location, will be overwritten 
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/reader2.py -l \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/spatial_cell_locations.txt\" -g \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/spatial_genes.txt\" -n \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/spatial_network.txt\" -e \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/expression_matrix.txt\" -o \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 10 2 2 -t 1e-10 -z none -s 100 -i 100"

# check and select hmrf
for(i in seq(10, 14, by = 2)) {
  viewHMRFresults2D(gobject = starmap_mini,
                    HMRFoutput = HMRF_spatial_genes,
                    k = 6, betas_to_view = i,
                    point_size = 2)
}
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 10"
#> Warning in system(command = result_command, intern = T): running command '/opt/
#> conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r
#> "../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore" -a
#> test -k 6 -b 10' had status 1
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 12"
#> Warning in system(command = result_command, intern = T): running command '/opt/
#> conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r
#> "../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore" -a
#> test -k 6 -b 12' had status 1

starmap_mini = addHMRF(gobject = starmap_mini,
                  HMRFoutput = HMRF_spatial_genes,
                  k = 6, betas_to_add = c(12),
                  hmrf_name = 'HMRF')
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 12"
#> Warning in system(command = result_command, intern = T): running command '/opt/
#> conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r
#> "../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore" -a
#> test -k 6 -b 12' had status 1

# visualize selected hmrf result
giotto_colors = Giotto:::getDistinctColors(6)
names(giotto_colors) = 1:6
spatPlot(gobject = starmap_mini, cell_color = 'HMRF_k6_b.12',
         point_size = 3, coord_fix_ratio = 1, cell_color_code = giotto_colors)

set.seed(seed = 2841)
cell_proximities = cellProximityEnrichment(gobject = starmap_mini,
                                           cluster_column = 'cell_types',
                                           spatial_network_name = 'Delaunay_network',
                                           adjust_method = 'fdr',
                                           number_of_simulations = 1000)
# barplot
cellProximityBarplot(gobject = starmap_mini, 
                     CPscore = cell_proximities, 
                     min_orig_ints = 5, min_sim_ints = 5)

## heatmap
cellProximityHeatmap(gobject = starmap_mini, CPscore = cell_proximities, 
                     order_cell_types = T, scale = T,
                     color_breaks = c(-1.5, 0, 1.5), 
                     color_names = c('blue', 'white', 'red'))

# network
cellProximityNetwork(gobject = starmap_mini, CPscore = cell_proximities, 
                     remove_self_edges = T, only_show_enrichment_edges = T)

# network with self-edges
cellProximityNetwork(gobject = starmap_mini, CPscore = cell_proximities,
                     remove_self_edges = F, self_loop_strength = 0.3,
                     only_show_enrichment_edges = F,
                     rescale_edge_weights = T,
                     node_size = 8,
                     edge_weight_range_depletion = c(1, 2),
                     edge_weight_range_enrichment = c(2,5))

# create cross section
starmap_mini = createCrossSection(starmap_mini,
                                method="equation",
                                equation=c(0,1,0,600),
                                extend_ratio = 0.6)

# show cross section
insertCrossSectionSpatPlot3D(starmap_mini, cell_color = 'leiden_clus',
                             axis_scale = 'cube',
                             point_size = 2)
#> create 3D plot

insertCrossSectionGenePlot3D(starmap_mini, expression_values = 'scaled',
                             axis_scale = "cube",
                             genes = "Slc17a7")

# for cell annotation
crossSectionPlot(starmap_mini,
                 point_size = 2, point_shape = "border",
                 cell_color = "leiden_clus")

crossSectionPlot3D(starmap_mini,
                   point_size = 2, cell_color = "leiden_clus",
                   axis_scale = "cube")
#> create 3D plot

# for gene expression
crossSectionGenePlot(starmap_mini,
                     genes = "Slc17a7",
                     point_size = 2,
                     point_shape = "border",
                     cow_n_col = 1.5,
                     expression_values = 'scaled')

crossSectionGenePlot3D(starmap_mini,
                       point_size = 2,
                       genes = c("Slc17a7"),
                       expression_values = 'scaled')

viewer_folder = paste0(temp_dir, '/', 'Mouse_cortex_viewer')

# select annotations, reductions and expression values to view in Giotto Viewer
exportGiottoViewer(gobject = starmap_mini, output_directory = viewer_folder,
                   factor_annotations = c('cell_types',
                                          'leiden_clus',
                                          'HMRF_k6_b.12'),
                   numeric_annotations = 'total_expr',
                   dim_reductions = c('umap'),
                   dim_reduction_names = c('umap'),
                   expression_values = 'scaled',
                   expression_rounding = 3,
                   overwrite_dir = T)
#> output directory already exits, files will be overwritten 
#> write cell and gene IDs
#> write physical centroid locations
#> write annotation data for:  cell_types 
#> write annotation data for:  leiden_clus 
#> write annotation data for:  HMRF_k6_b.12 
#> write annotation data for:  total_expr 
#> write annotation data for:  umap  for  umap 
#> write expression values 
#> finished writing giotto viewer files to ../results//Mouse_cortex_viewer 
#> 
#> ================================================================ 
#> Next steps. Please manually run the following in a SHELL terminal: 
#> ================================================================ 
#> cd  ../results//Mouse_cortex_viewer 
#> giotto_setup_image --require-stitch=n --image=n --image-multi-channel=n --segmentation=n --multi-fov=n --output-json=step1.json 
#> smfish_step1_setup -c step1.json 
#> giotto_setup_viewer --num-panel=2 --input-preprocess-json=step1.json --panel-1=PanelPhysicalSimple --panel-2=PanelTsne --output-json=step2.json --input-annotation-list=annotation_list.txt 
#> smfish_read_config -c step2.json -o test.dec6.js -p test.dec6.html -q test.dec6.css 
#> giotto_copy_js_css --output . 
#> python3 -m http.server 
#> ================================================================ 
#> 
#> Finally, open your browser, navigate to http://localhost:8000/. Then click on the file test.dec6.html to see the viewer. 
#> 
#> 
#> For more information, http://spatialgiotto.rc.fas.harvard.edu/giotto.viewer.setup3.html